This script is to compare model fits between 1) the old model (unpooled K), 2) the pooled-K model we developed, with pooled and logged K, and 3) a non-centered, pooled-K model.
Summary: the new pooled-logK model is still really struggling for two sites: BIGC and CARI. ChatGPT and other resources suggested a non-centered model structure to avoid ‘funnelling’ for sigma and mu: instead of drawing each logU or logK directly from a normal distribution with unknown σ, the non-centered version draws a standard-normal z first, then scales and shifts it by μ and σ.
For BIGC, the new data is part of the problem: the earlier, unpooled-K model is also struggling. The non-centered model works well. Solution: find days that might be causing the problems, remove from the dataset, try all models again.
For CARI, the original unpooled-K model runs well, but not the initial pooled-K model. The uncentered pooled-K model works well.
Parameter values are similar across all three models, which is a good sign.
For the other four sites, CUPE, PRIN, SYCA, WLOU:
The original pooled-K model does very well for CUPE and SYCA, and decently well for PRIN and WLOU. More iterations will probably bring all r-hats to 1.
The non-centered pooled-K model is supposed to improve model performance across the board, but it does not. As noted, it worked really well for BIGC and CARI, decently well but more slowly (vs the other pooled-K model) for most of the others (CUPE, PRIN, SYCA), and poorly for WLOU. Again, parameter estimates are similar across all three models (unpooled, pooled logK, uncentered pooled logK), so not much to choose from there.
BIGC probably has problem days, removing them will hopefully improve the modeling. CARI does well w. the non-centered model and struggles with our initial pooled-K (‘centered’) model. WLOU does well w. the ‘centered’ model and struggles with the non-centered model. Thoughts on how to proceed?
########## Load model data _____________________________________________________# The same model data should work for both pooled and unpooled modelsbigc.modeldata <-readRDS(here("N_uptake_NEON/data/model_datalist/bigc.data.rds"))cari.modeldata <-readRDS(here("N_uptake_NEON/data/model_datalist/cari.data.rds"))cupe.modeldata <-readRDS(here("N_uptake_NEON/data/model_datalist/cupe.data.rds"))prin.modeldata <-readRDS(here("N_uptake_NEON/data/model_datalist/prin.data.rds"))syca.modeldata <-readRDS(here("N_uptake_NEON/data/model_datalist/syca.data.rds"))wlou.modeldata <-readRDS(here("N_uptake_NEON/data/model_datalist/wlou.data.rds"))
Big Creek (BIGC)
Unpooled K
Model finished in 20 min 10 sec BIGC currently struggles with both the old and new model, and clearly needs some ‘problem’ days removed.
Show the code
########## Load unpooled model fit ______________________________________________________bigc.fit_unpooledK <-readRDS(here("N_uptake_NEON/data/model_fits/unpooledK/bigc.fit.rds"))
Inference for Stan model: anon_model.
4 chains, each with iter=1000; warmup=500; thin=1;
post-warmup draws per chain=500, total post-warmup draws=2000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
sigma 0.06 0.00 0.00 0.06 0.06 0.06 0.06 0.06 2049 1.00
sigma_U 0.07 0.01 0.03 0.02 0.05 0.07 0.09 0.15 15 1.28
b0 -7.30 0.13 0.94 -9.29 -7.92 -7.25 -6.65 -5.55 53 1.07
b1 0.41 0.01 0.08 0.26 0.36 0.41 0.47 0.59 53 1.07
Samples were drawn using NUTS(diag_e) at Tue Oct 28 20:23:46 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).
Examine MCMC chains
Show the code
posterior <-as.array(bigc.fit_unpooledK) # from rstan packagecolor_scheme_set("viridis")mcmc_trace(posterior, pars =c("sigma", "sigma_U", "b0", "b1"))
Show the code
# highlights one chain (uses points not lines)mcmc_trace_highlight(posterior, pars =c("sigma", "sigma_U", "b0", "b1"), highlight =2, size =2)
Pooled logK
Model finished in 21 min 54 sec
Show the code
########## Load pooled logK model fit ______________________________________________________bigc.fit_pooledK.log <-readRDS(here("N_uptake_NEON/data/model_fits/pooledK_logK/bigc.fit.rds"))
Print pooled logK parameter fit
Show the code
########## Load pooled logK model fit ______________________________________________________print(bigc.fit_pooledK.log, pars=c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"))
Inference for Stan model: anon_model.
4 chains, each with iter=1000; warmup=500; thin=1;
post-warmup draws per chain=500, total post-warmup draws=2000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
sigma 0.06 0.00 0.00 0.06 0.06 0.06 0.06 0.06 746 1.00
sigma_U 0.08 0.01 0.03 0.02 0.05 0.07 0.10 0.14 15 1.30
sigma_logK 0.44 0.01 0.07 0.30 0.39 0.44 0.49 0.59 74 1.04
logK_mean 0.83 0.01 0.09 0.65 0.78 0.83 0.89 1.00 75 1.04
b0 -7.39 0.36 1.29 -10.45 -8.06 -7.20 -6.45 -5.33 13 1.32
b1 0.42 0.03 0.11 0.23 0.34 0.40 0.48 0.69 13 1.32
Samples were drawn using NUTS(diag_e) at Tue Oct 28 09:53:24 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).
Examine MCMC chains
Show the code
posterior <-as.array(bigc.fit_pooledK.log) # from rstan packagecolor_scheme_set("viridis")mcmc_trace(posterior, pars =c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"))
Show the code
# highlights one chain (uses points not lines)mcmc_trace_highlight(posterior, pars =c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"), highlight =4, size =2)
Show the code
##### Examine distributions of struggling params + # # gets fit into a dataframe - easier for the bayesplotting# posteriorMA <- as_draws_matrix(posterior)# sigma <- posteriorMA[, "sigma"]# # # Find all conc_tilde columns# conc_cols <- grep("^conc_tilde\\[", colnames(posteriorMA), value = TRUE)# # cors <- sapply(conc_cols, function(nm) {# x <- posteriorMA[, nm]# if (is.numeric(x) && all(!is.na(x))) {# cor(x, sigma)# } else {# NA# }# })# # cors <- cors[!is.na(cors)]# top <- names(sort(abs(cors), decreasing = TRUE))[1]# # cat("Most correlated:", top, " (r =", round(cors[top], 3), ")\n")# # # bayesplot::mcmc_scatter(as.array(bigc.fit_pooledK.log),# pars = c(top, "sigma"),# transform = list(sigma = "log"), # alpha = 0.3)# # # # scatter_sigma_pool <- mcmc_scatter(# posterior, # pars = c("sigma", "conc_tilde[10,20]"), # transform = list(sigma = "log"), # size = 1# )# # scatter_sigma_pool#
Uncentered and pooled
A second, non-centered pooled model worked better: - Model finished in 7 min 45 sec
Show the code
########## Load pooled logK model fit ______________________________________________________bigc.fit2_pooledK.log <-readRDS(here("N_uptake_NEON/data/model_fits/pooledK_logK/bigc.fit2.rds"))
Print non-centered pooled logK parameter fit
Show the code
########## Load pooled logK model fit ______________________________________________________print(bigc.fit2_pooledK.log, pars=c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"))
Inference for Stan model: anon_model.
4 chains, each with iter=1000; warmup=500; thin=1;
post-warmup draws per chain=500, total post-warmup draws=2000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
sigma 0.06 0.00 0.00 0.06 0.06 0.06 0.06 0.06 2320 1.00
sigma_U 0.04 0.00 0.03 0.00 0.02 0.04 0.06 0.12 765 1.00
sigma_logK 0.44 0.00 0.08 0.30 0.38 0.43 0.49 0.59 266 1.02
logK_mean 0.82 0.00 0.09 0.64 0.76 0.82 0.88 0.98 399 1.01
b0 -6.99 0.02 0.97 -8.85 -7.61 -6.99 -6.32 -5.16 1682 1.00
b1 0.38 0.00 0.09 0.22 0.32 0.38 0.44 0.55 1674 1.00
Samples were drawn using NUTS(diag_e) at Tue Oct 28 23:35:13 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).
Examine MCMC chains
Show the code
posterior <-as.array(bigc.fit2_pooledK.log) # from rstan packagecolor_scheme_set("viridis")mcmc_trace(posterior, pars =c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"))
Show the code
# highlights one chain (uses points not lines) - all chains converged well# mcmc_trace_highlight(posterior, pars = c("sigma", "sigma_U", "b0", "b1"), highlight = 4, size = 2)
Caribou Creek (CARI)
Unpooled K
Model finished in 3 min 36 sec
Show the code
########## Load unpooled model fit ______________________________________________________cari.fit_unpooledK <-readRDS(here("N_uptake_NEON/data/model_fits/unpooledK/cari.fit.rds"))
Inference for Stan model: anon_model.
4 chains, each with iter=1000; warmup=500; thin=1;
post-warmup draws per chain=500, total post-warmup draws=2000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
sigma 0.13 0.00 0.00 0.12 0.13 0.13 0.13 0.13 4363 1
sigma_U 0.40 0.00 0.03 0.34 0.37 0.40 0.42 0.47 579 1
b0 -1.95 0.02 0.53 -3.00 -2.27 -1.94 -1.59 -0.94 1122 1
b1 0.12 0.00 0.05 0.03 0.09 0.12 0.15 0.21 1111 1
Samples were drawn using NUTS(diag_e) at Tue Oct 28 20:37:16 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).
Examine MCMC chains
Show the code
posterior <-as.array(cari.fit_unpooledK) # from rstan packagecolor_scheme_set("viridis")mcmc_trace(posterior, pars =c("sigma", "sigma_U", "b0", "b1"))
Show the code
# highlights one chain (uses points not lines) - not needed when chains look good# mcmc_trace_highlight(posterior, pars = c("sigma", "sigma_U", "b0", "b1"), highlight = 2, size = 2)
Pooled logK
Model finished in 5 min 15 sec
Show the code
########## Load pooled logK model fit ______________________________________________________cari.fit_pooledK.log <-readRDS(here("N_uptake_NEON/data/model_fits/pooledK_logK/cari.fit.rds"))
A second, non-centered pooled model composition worked better:
Show the code
########## Load pooled logK model fit ______________________________________________________cari.fit2_pooledK.log <-readRDS(here("N_uptake_NEON/data/model_fits/pooledK_logK/cari.fit2.rds"))
Print non-centered pooled logK parameter fit
Show the code
########## Load pooled logK model fit ______________________________________________________print(cari.fit2_pooledK.log, pars=c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"))
Inference for Stan model: anon_model.
4 chains, each with iter=1000; warmup=500; thin=1;
post-warmup draws per chain=500, total post-warmup draws=2000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
sigma 0.12 0.00 0.00 0.12 0.12 0.12 0.12 0.13 2890 1
sigma_U 0.38 0.00 0.03 0.32 0.36 0.38 0.40 0.45 926 1
sigma_logK 0.06 0.00 0.04 0.00 0.02 0.05 0.08 0.17 780 1
logK_mean -1.00 0.00 0.05 -1.10 -1.04 -1.00 -0.97 -0.91 1202 1
b0 -2.32 0.01 0.50 -3.32 -2.65 -2.32 -1.97 -1.36 1201 1
b1 0.15 0.00 0.04 0.07 0.12 0.15 0.18 0.24 1184 1
Samples were drawn using NUTS(diag_e) at Tue Oct 28 23:52:03 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).
Examine MCMC chains
Show the code
posterior <-as.array(cari.fit2_pooledK.log) # from rstan packagecolor_scheme_set("viridis")mcmc_trace(posterior, pars =c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"))
Show the code
# highlights one chain (uses points not lines) - all chains converged well# mcmc_trace_highlight(posterior, pars = c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"), highlight = 4, size = 2)
Rio Cupeyes (CUPE)
Unpooled K
Model finished in 9 min 16 sec
Show the code
########## Load unpooled model fit ______________________________________________________cupe.fit_unpooledK <-readRDS(here("N_uptake_NEON/data/model_fits/unpooledK/cupe.fit.rds"))
Inference for Stan model: anon_model.
4 chains, each with iter=1000; warmup=500; thin=1;
post-warmup draws per chain=500, total post-warmup draws=2000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
sigma 0.08 0.00 0.00 0.08 0.08 0.08 0.08 0.08 2730 1
sigma_U 0.30 0.00 0.01 0.27 0.29 0.30 0.30 0.32 2573 1
b0 -3.08 0.01 0.62 -4.29 -3.50 -3.09 -2.67 -1.83 1816 1
b1 0.13 0.00 0.06 0.02 0.09 0.13 0.17 0.24 1821 1
Samples were drawn using NUTS(diag_e) at Tue Oct 28 20:51:59 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).
Examine MCMC chains
Show the code
posterior <-as.array(cupe.fit_unpooledK) # from rstan packagecolor_scheme_set("viridis")mcmc_trace(posterior, pars =c("sigma", "sigma_U", "b0", "b1"))
Show the code
# highlights one chain (uses points not lines) - good convergence, did not need# mcmc_trace_highlight(posterior, pars = c("sigma", "sigma_U", "b0", "b1"), highlight = 2, size = 2)
Pooled logK
Model finished in 18 min 56 sec
Show the code
########## Load pooled logK model fit ______________________________________________________cupe.fit_pooledK.log <-readRDS(here("N_uptake_NEON/data/model_fits/pooledK_logK/cupe.fit.rds"))
A second, non-centered pooled model composition also worked for CUPE, but not quite as well as the original pooled-K - and took a VERY long time, maybe partly due to other draws on computation power/ bandwidth
Show the code
########## Load pooled logK model fit ______________________________________________________cupe.fit2_pooledK.log <-readRDS(here("N_uptake_NEON/data/model_fits/pooledK_logK/cupe.fit2.rds"))
Print non-centered pooled logK parameter fit
Show the code
########## Load pooled logK model fit ______________________________________________________print(cupe.fit2_pooledK.log, pars=c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"))
Inference for Stan model: anon_model.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
sigma 0.08 0.00 0.00 0.08 0.08 0.08 0.08 0.08 5206 1.00
sigma_U 0.29 0.00 0.01 0.27 0.28 0.29 0.29 0.31 1015 1.00
sigma_logK 0.21 0.00 0.01 0.18 0.20 0.21 0.22 0.23 1782 1.00
logK_mean 1.40 0.00 0.01 1.37 1.39 1.40 1.40 1.42 2808 1.00
b0 -3.10 0.03 0.61 -4.28 -3.50 -3.10 -2.69 -1.94 466 1.02
b1 0.13 0.00 0.06 0.03 0.09 0.13 0.17 0.24 458 1.02
Samples were drawn using NUTS(diag_e) at Wed Oct 29 01:16:26 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).
Examine MCMC chains
Show the code
posterior <-as.array(cupe.fit2_pooledK.log) # from rstan packagecolor_scheme_set("viridis")mcmc_trace(posterior, pars =c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"))
Show the code
# highlights one chain (uses points not lines) - all chains converged well# mcmc_trace_highlight(posterior, pars = c("sigma", "sigma_U", "b0", "b1"), highlight = 4, size = 2)
Pringle Creek (PRIN)
Unpooled K
Model finished in 3 min 16 sec
Show the code
########## Load unpooled model fit ______________________________________________________prin.fit_unpooledK <-readRDS(here("N_uptake_NEON/data/model_fits/unpooledK/prin.fit.rds"))
Inference for Stan model: anon_model.
4 chains, each with iter=1000; warmup=500; thin=1;
post-warmup draws per chain=500, total post-warmup draws=2000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
sigma 0.29 0.00 0.00 0.28 0.29 0.29 0.29 0.30 2841 1
sigma_U 0.49 0.00 0.04 0.42 0.46 0.49 0.51 0.57 976 1
b0 -11.87 0.05 1.69 -15.17 -13.01 -11.83 -10.69 -8.58 980 1
b1 0.98 0.00 0.15 0.68 0.87 0.97 1.08 1.27 988 1
Samples were drawn using NUTS(diag_e) at Tue Oct 28 21:02:43 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).
Examine MCMC chains
Show the code
posterior <-as.array(prin.fit_unpooledK) # from rstan packagecolor_scheme_set("viridis")mcmc_trace(posterior, pars =c("sigma", "sigma_U", "b0", "b1"))
Show the code
# highlights one chain (uses points not lines) - chains converged, did not need# mcmc_trace_highlight(posterior, pars = c("sigma", "sigma_U", "b0", "b1"), highlight = 2, size = 2)
Pooled logK
Model finished in 3 min 40 sec
Show the code
########## Load pooled logK model fit ______________________________________________________prin.fit_pooledK.log <-readRDS(here("N_uptake_NEON/data/model_fits/pooledK_logK/prin.fit.rds"))
A second, non-centered pooled model composition also worked, but took longer to converge:
Show the code
########## Load pooled logK model fit ______________________________________________________prin.fit2_pooledK.log <-readRDS(here("N_uptake_NEON/data/model_fits/pooledK_logK/prin.fit2.rds"))
Print non-centered pooled logK parameter fit
Show the code
########## Load pooled logK model fit ______________________________________________________print(prin.fit2_pooledK.log, pars=c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"))
Inference for Stan model: anon_model.
4 chains, each with iter=1000; warmup=500; thin=1;
post-warmup draws per chain=500, total post-warmup draws=2000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
sigma 0.28 0.00 0.00 0.28 0.28 0.28 0.29 0.29 2722 1.00
sigma_U 0.47 0.00 0.04 0.40 0.44 0.47 0.50 0.55 609 1.00
sigma_logK 0.39 0.01 0.13 0.12 0.31 0.39 0.48 0.66 152 1.02
logK_mean -0.74 0.01 0.11 -0.98 -0.81 -0.73 -0.66 -0.56 182 1.01
b0 -11.32 0.08 1.66 -14.54 -12.45 -11.33 -10.17 -8.09 434 1.00
b1 0.93 0.01 0.15 0.64 0.82 0.93 1.03 1.21 425 1.00
Samples were drawn using NUTS(diag_e) at Wed Oct 29 01:30:56 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).
Examine MCMC chains
Show the code
posterior <-as.array(prin.fit2_pooledK.log) # from rstan packagecolor_scheme_set("viridis")mcmc_trace(posterior, pars =c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"))
Show the code
# highlights one chain (uses points not lines) - all chains converged fairly well with a couple of wobbles in the logK hyperparametersmcmc_trace_highlight(posterior, pars =c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"), highlight =1, size =2, alpha=0.5)
########## Load unpooled model fit ______________________________________________________syca.fit_unpooledK <-readRDS(here("N_uptake_NEON/data/model_fits/unpooledK/syca.fit.rds"))
Inference for Stan model: anon_model.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
sigma 0.43 0.00 0.01 0.42 0.43 0.43 0.44 0.45 5439 1
sigma_U 0.75 0.00 0.07 0.63 0.70 0.74 0.79 0.89 2799 1
b0 32.54 0.04 2.91 26.78 30.60 32.60 34.50 38.27 4423 1
b1 -3.77 0.01 0.33 -4.43 -4.00 -3.78 -3.55 -3.12 4358 1
Samples were drawn using NUTS(diag_e) at Tue Oct 28 21:20:27 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).
Examine MCMC chains
Show the code
posterior <-as.array(syca.fit_unpooledK) # from rstan packagecolor_scheme_set("viridis")mcmc_trace(posterior, pars =c("sigma", "sigma_U", "b0", "b1"))
Show the code
# highlights one chain (uses points not lines) - chains converged, did not need# mcmc_trace_highlight(posterior, pars = c("sigma", "sigma_U", "b0", "b1"), highlight = 2, size = 2)
Pooled logK
Model finished in 2 min 22 sec
Show the code
########## Load pooled logK model fit ______________________________________________________syca.fit_pooledK.log <-readRDS(here("N_uptake_NEON/data/model_fits/pooledK_logK/syca.fit.rds"))
Inference for Stan model: anon_model.
4 chains, each with iter=1000; warmup=500; thin=1;
post-warmup draws per chain=500, total post-warmup draws=2000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
sigma 0.42 0.00 0.01 0.41 0.42 0.42 0.43 0.43 3119 1.00
sigma_U 0.62 0.00 0.05 0.53 0.58 0.62 0.66 0.73 2148 1.00
sigma_logK 0.18 0.00 0.03 0.12 0.16 0.18 0.20 0.25 465 1.01
logK_mean 2.08 0.00 0.03 2.02 2.06 2.08 2.11 2.15 949 1.00
b0 26.15 0.05 2.44 21.30 24.51 26.17 27.75 30.90 2471 1.00
b1 -3.00 0.01 0.28 -3.54 -3.18 -3.01 -2.81 -2.45 2440 1.00
Samples were drawn using NUTS(diag_e) at Tue Oct 28 11:23:26 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).
Examine MCMC chains
Show the code
posterior <-as.array(syca.fit_pooledK.log) # from rstan packagecolor_scheme_set("viridis")mcmc_trace(posterior, pars =c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"))
Show the code
# highlights one chain (uses points not lines) - chains converged, did not need# mcmc_trace_highlight(posterior, pars = c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"), highlight = 4, size = 2)
Uncentered and pooled
Model finished in 8 min 42 sec
A second, non-centered pooled model composition also worked, but took longer to converge:
Show the code
########## Load pooled logK model fit ______________________________________________________syca.fit2_pooledK.log <-readRDS(here("N_uptake_NEON/data/model_fits/pooledK_logK/syca.fit2.rds"))
Print non-centered pooled logK parameter fit
Show the code
########## Load pooled logK model fit ______________________________________________________print(syca.fit2_pooledK.log, pars=c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"))
Inference for Stan model: anon_model.
4 chains, each with iter=1000; warmup=500; thin=1;
post-warmup draws per chain=500, total post-warmup draws=2000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
sigma 0.42 0.00 0.01 0.41 0.42 0.42 0.43 0.43 4186 1.00
sigma_U 0.62 0.00 0.05 0.53 0.59 0.62 0.66 0.73 571 1.01
sigma_logK 0.18 0.00 0.03 0.12 0.16 0.18 0.20 0.25 884 1.00
logK_mean 2.08 0.00 0.03 2.02 2.06 2.08 2.11 2.15 942 1.00
b0 25.89 0.14 2.43 21.16 24.22 25.95 27.45 30.55 298 1.00
b1 -2.97 0.02 0.28 -3.51 -3.15 -2.98 -2.78 -2.44 311 1.00
Samples were drawn using NUTS(diag_e) at Wed Oct 29 01:57:05 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).
Examine MCMC chains
Show the code
posterior <-as.array(syca.fit2_pooledK.log) # from rstan packagecolor_scheme_set("viridis")mcmc_trace(posterior, pars =c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"))
Show the code
# highlights one chain (uses points not lines) - mb a little wobbly in the betas(?)mcmc_trace_highlight(posterior, pars =c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"), highlight =1, size =2, alpha=0.5)
########## Load unpooled model fit ______________________________________________________wlou.fit_unpooledK <-readRDS(here("N_uptake_NEON/data/model_fits/unpooledK/wlou.fit.rds"))
Inference for Stan model: anon_model.
4 chains, each with iter=4000; warmup=2000; thin=1;
post-warmup draws per chain=2000, total post-warmup draws=8000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
sigma 0.05 0.00 0.00 0.05 0.05 0.05 0.05 0.05 9714 1
sigma_U 0.64 0.00 0.04 0.57 0.62 0.64 0.67 0.72 4044 1
b0 -6.74 0.01 0.89 -8.48 -7.34 -6.74 -6.14 -4.98 4338 1
b1 0.43 0.00 0.09 0.25 0.37 0.43 0.49 0.61 4341 1
Samples were drawn using NUTS(diag_e) at Tue Oct 28 21:41:19 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).
Examine MCMC chains
Show the code
posterior <-as.array(wlou.fit_unpooledK) # from rstan packagecolor_scheme_set("viridis")mcmc_trace(posterior, pars =c("sigma", "sigma_U", "b0", "b1"))
Show the code
# highlights one chain (uses points not lines) - chains converged well, did not need# mcmc_trace_highlight(posterior, pars = c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"), highlight = 2, size = 2)
Pooled logK
Model finished in 6 min 27 sec
Show the code
########## Load pooled logK model fit ______________________________________________________wlou.fit_pooledK.log <-readRDS(here("N_uptake_NEON/data/model_fits/pooledK_logK/wlou.fit.rds"))
A second, non-centered pooled model composition did not work well for WLOU:
Show the code
########## Load pooled logK model fit ______________________________________________________wlou.fit2_pooledK.log <-readRDS(here("N_uptake_NEON/data/model_fits/pooledK_logK/wlou.fit2.rds"))
Print non-centered pooled logK parameter fit
Show the code
########## Load pooled logK model fit ______________________________________________________print(wlou.fit2_pooledK.log, pars=c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"))
Inference for Stan model: anon_model.
4 chains, each with iter=2500; warmup=1250; thin=1;
post-warmup draws per chain=1250, total post-warmup draws=5000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
sigma 0.05 0.00 0.00 0.05 0.05 0.05 0.05 0.05 3066 1.00
sigma_U 0.64 0.00 0.04 0.57 0.61 0.64 0.66 0.71 1549 1.00
sigma_logK 0.48 0.00 0.07 0.35 0.43 0.48 0.53 0.65 330 1.01
logK_mean 1.14 0.00 0.08 0.95 1.09 1.15 1.19 1.27 395 1.01
b0 -6.84 0.03 0.89 -8.63 -7.43 -6.84 -6.23 -5.16 964 1.00
b1 0.44 0.00 0.09 0.27 0.38 0.44 0.51 0.63 969 1.00
Samples were drawn using NUTS(diag_e) at Mon Nov 24 13:50:48 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).
Examine MCMC chains
Show the code
posterior <-as.array(wlou.fit2_pooledK.log) # from rstan packagecolor_scheme_set("viridis")mcmc_trace(posterior, pars =c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"))
Show the code
# highlights one chain (uses points not lines) mcmc_trace_highlight(posterior, pars =c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"), highlight =1, size =2, alpha=0.5)